# -------------------------------------------------------------------
# Purpose: Creates Figure B9
# Author:  Max Posch, 25/07/2025
# Usage:   Source this script to generate all panels of the figure
# -------------------------------------------------------------------
# Check that required paths exist
stopifnot(dir.exists(pdataraw))
stopifnot(dir.exists(pdataanalysis))
stopifnot(dir.exists(poutputappendix))

p_unload(tidytable)


## Load data
load(file.path(pdataanalysis, "countyLevel19001940.RData"))
load(file.path(pdataraw, "US_tl2000_cnty_shp_1790_2000.Rdata"))
load(file.path(pdataraw, "US_tl2000_state_shp_1790_2000.Rdata"))


## to address warning about old-style crs
st_crs(shp2000[["1900"]]) <- st_crs(shp2000[["1900"]])
st_crs(state_shp2000[["1900"]]) <- st_crs(state_shp2000[["1900"]])


## simplify shapes
cnty1900 <- shp2000[["1900"]] %>%
    filter(ICPSRSTI != 82, ICPSRSTI != 81, ICPSRSTI != 0) %>%
    select(gisjoin_1900 = GISJOIN, geometry) %>%
    rmapshaper::ms_simplify(keep = 0.05, keep_shapes = TRUE)

state1900 <- state_shp2000[["1900"]] %>%
        filter(ICPSRST != 82, ICPSRST != 81) %>%
        select(geometry) %>%
        rmapshaper::ms_simplify(keep = 0.05, keep_shapes = TRUE)


## Subset data
d <- na.omit(countyLevel19001940[, .(gisjoin_1900, iv_lo_entropy_namelast_mp_adjp_fe_immig_w, iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_w, statefip, year)])


## residualize
d[, iv_lo_entropy_namelast_mp_adjp_fe_immig_r_w := as.numeric(resid(feols(iv_lo_entropy_namelast_mp_adjp_fe_immig_w ~ iv_lo_sum_n_namelast_mp_adjp_fe_immig_tr_w | gisjoin_1900 + statefip^year, d)))]


## Create plots for multiple years
years_to_plot <- c(1900, 1905, 1910, 1915, 1920, 1925, 1930, 1940)
plot_files <- sapply(years_to_plot, function(yr) {
    cat(paste0("Creating plot for year ", yr, "...\n"))
    plotname <- create_figureB09(yr)
    if(!is.null(plotname)) {
        cat("Figure saved to:", plotname, "\n")
    }
    return(plotname)
})
p_load(tidytable)
